Dynamics of the Wang-Landau algorithm and complexity of rare events 
for the three-dimensional bimodal Ising spin glass 
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We investigate the performance of flat-histogram methods based on a multicanonical ensemble 
and the Wang-Landau algorithm for the three-dimensional ±J spin glass by measuring round-trip 
times in the energy range between the zero-temperature ground state and the state of highest energy. 
Strong sample-to-sample variations are found for fixed system size and the distribution of round-trip 
times follows a fat-tailed Frechet extremal value distribution. Rare events in the fat tails of these 
distributions corresponding to extremely slowly equilibrating spin glass realizations dominate the 
calculations of statistical averages. While the typical round-trip time scales exponential as expected 
for this NP-hard problem, we find that the average round-trip time is no longer well-defined for 
systems with N >8 3 spins. We relate the round-trip times for multicanonical sampling to intrinsic 
properties of the energy landscape and compare with the numerical effort needed by the genetic 
Cluster-Exact Approximation to calculate the exact ground state energies. For systems with N > 8 3 
spins the simulation of these rare events becomes increasingly hard. For N > 14 3 there are samples 
where the Wang-Landau algorithm fails to find the true ground state within reasonable simulation 
times. We expect similar behavior for other algorithms based on multicanonical sampling. 



PACS numbers: 02.70.Rr, 75.10.Hk, 64.60.Cn 



I. INTRODUCTION 

The three-dimensional ± J Ising spin glass has been ex- 
tensively studied^ as a prototype system which exhibits 
a finite temperature second order phase transition to a 
slowly equilibrating glassy phase The simulation of 
such a system with conventional Monte Carlo methods is 
slowed down by long relaxation times in the spin-glass 
phase. This problem has been addressed by a num- 
ber of algorithmic developments such as the multicanon- 
ical method^ simulated and parallel tempering, 5 broad 
histograms^ and transition matrix Monte Carlo^In order 
to speed up equilibration most of these methods aim at 
broadening the energy range sampled within the Monte 
Carlo (MC) simulations from the sharply peaked distri- 
bution of canonical sampling at a fixed temperature. 

Recently, Wang and Landau introduced a new algo- 
rithm which systematically calculates an estimate of the 
density of states and iteratively converges to sampling a 
flat histogram in energy^ The Wang-Landau (WL) al- 
gorithm simulates a biased random walk in configura- 
tion space. The bias depends only on the total energy 
of a configuration and is defined by a statistical ensem- 
ble with weights w(E). For this ensemble the transition 
probabilities are given by the Metropolis scheme 



p(Ei -> E 2 ) = min ( W ^?\ , 1 
1 w[Ei) 



(1) 



The equilibrium distribution of the energy in this en- 
semble is n w (E) oc w(E)g(E) where g(E) is the density 
of states. By setting the weights w(E) oc l/g(E) the 
WL algorithm aims at sampling a flat histogram in en- 
ergy. The crucial feature of the WL algorithm however is 



that the simulated ensemble is dynamically modified dur- 
ing the course of the simulation: after every spin update 
the current estimate of the density of states is multiplied 
by a modification factor / and the ensemble weights are 
analogously updated. The modification factor is itera- 
tively reduced to 1 whenever the sampled energy his- 
togram is close to the expected equilibrium distribution, 
that is when the histogram is "flat" within a given range. 
The dynamic modification of the ensemble allows to push 
the random walker towards the low entropy states in the 
initial stages of the algorithm while ensuring that it con- 
verges to a flat-histogram/multicanonical ensemble in the 
final stages of the computation. 

In this paper, we determine the performance of the 
Wang-Landau algorithm for the three-dimensional ±J 
Ising spin glass for both stages of the algorithm. First, 
we study the dynamic behavior of the algorithm in the 
initial stages by considering its ability to find the ground- 
state energy of a number of three-dimensional spin glass 
samples which is a well known NP-hard problem £ We 
compare the obtained ground state energies to exact re- 
sults calculated with the genetic Cluster-Exact approxi- 
mation (CEA ) 14115 . While the WL algorithm reproduces 
the exact ground-state energy for small systems, we find 
that for moderately large systems (N > 14 3 ) the WL al- 
gorithm does not find the exact ground-state energy for 
few spin-glass samples. Even when restricting the simu- 
lated energy bin around the known ground-state energy 
the algorithm does not find the lowest energy state within 
a reasonable number of sweeps (A swe eps ~ 10 7 ). The ge- 
netic CEA gives a superior performance to find ground 
state energies, but does not give the full thermodynamic 
information. Second, we investigate the asymptotic be- 
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havior of the WL algorithm by measuring round-trip 
times in energy for the converged ensemble. The round- 
trip time gives a direct estimate of the equilibration time 
for the multicanonical ensemble. The asymptotic scal- 
ing of the Wang-Landau algorithm therefore also reflects 
the performance of other flat-histogram methods based 
on multicanonical sampling, such as the multicanonical 
method, 4 broad histograms 6 or transition matrix Monte 
Carlo.' We find large sample-to-sample variations of the 
round-trip times which can be described by fat-tailed 
Frechet extremal-value distributions. We discuss impor- 
tant implications for statistical sample averaging which 
are caused by the rare events in the fat tails of these 
distributions. The intrinsic character of the observed 
extremal-value distributions is demonstrated. Finally it 
is shown that these distributions scale exponentially with 
the linear system size. In comparison, we find that the 
computational effort of the genetic CEA is also correlated 
to the density of states, but is less sensitive to the low- 
energy landscape than the WL algorithm. Finally, we 
discuss these limitations of multicanonical sampling by 
measuring the local diffusivity of the random walker in 
energy. We find a pronounced minimum of the diffusivity 
near the ground-state energy which is symptomatic for 
the entropic barrier which slows down the equilibration 
of the random walker. 



II. DYNAMIC PERFORMANCE 

To study the dynamic performance of the Wang- 
Landau algorithm we have tested its ability to find the 
ground-state energy of a number of spin glass samples 
where the exact ground state is known. Finding the 
ground-state energy of the ±J spin-glass is an exten- 
sively studied problem, both by physicists and com- 
puter scientists. For the two-dimensional ± J spin glass 
the full density of states can be calculated with poly- 
nomial effort^Siii However, for the three-dimensional 
±J spin glass as well as the two-dimensional ±J spin 
glass with an external field, even the problem of finding 
the ground-state energy has been shown to be NP-hard. 9 
The direct calculation of the ground-state energy of spin 
glass samples using sophisticated exact branch-and-cut 
algorithms 12,13 is therefore limited to rather small sys- 
tems. In this study, we applied a combination of a genetic 
algorithm with Cluster-Exact Approximation (CEA) to 
determine ground state energies piii^ The algorithm, al- 
though based on heuristics, is able to find true ground 
states of samples up to size N = 16 3 = 4096 in reason- 
able, although exponentially growing time— 

We use the results obtained from this genetic + CEA 
approach to check the accuracy of the ground-state ener- 
gies found with the WL algorithm for samples up to size 
N = 14 3 . For system size, N < 6 3 , the density of states 
of each sample was calculated without energy binning 
by 50 independent runs. For all runs we find that the 
WL algorithm gives the exact ground-state energy. For 



N = 8 3 , for some few samples (17 out of 1000) the exact 
ground-state energy was found only by running exten- 
sive runs after the comparison with the results from the 
heuristic approach revealed that the true ground state 
energy has not been found. 

For larger systems with TV = 10 3 and more spins we 
have restricted the WL simulations to the smallest al- 
lowed energy range around the ground-state energy cal- 
culated by the genetic CEA. In order to keep ergodicity 
the energy bin has to be larger than 

AE/J > 4L d ~ 1 , (2) 

to assure that two domain walls can be inserted which 
thus enables moves which subsequently flip all spins and 
ergodicity within the energy bin is given. 

Each sample was simulated by six independent runs. 
For all samples with N = 10 3 spins the WL algorithm 
found true ground-state energies. However, for three 
samples not all runs gave the true ground-state energy, 
but sampled energies down to the first excited state only 
(within (2.8 ± 0.2) x 10 7 MC sweeps). For the samples 
with N = 12 3 spins we find similar results. For all sam- 
ples there is at least one run which finds the true ground- 
state energy within some (2.9 ± 0.1) x 10 7 MC sweeps. 
However, for 9 out of 10 samples the WL algorithm does 
not converge to sampling a flat histogram in the full en- 
ergy range. For the samples with N — 14 3 spins the 
WL algorithms finds the true ground-state energy for 4 
samples only. For 6 out of 10 samples the algorithm did 
not sample the exact ground-state energy once within 
(3.0 ± 0.1) x 10 7 MC sweeps for all independent runs. 
Furthermore, there is one sample where even the first 
excited state is not found within the given number of 
sweeps. For all samples the simulated ensemble does not 
converge towards the multicanonical ensemble sampling 
a flat energy histogram and the WL algorithm gets stuck. 

III. ASYMPTOTIC PERFORMANCE 

We now turn to the asymptotic behavior of the WL 
algorithm and flat-histogram methods in general which 
we determine by measuring the round-trip times in en- 
ergy of the simulated random walker for the converged 
flat-histogram ensemble. Here the round-trip time corre- 
sponds to the number of single-spin flips needed to get 
from a configuration with the ground-state energy to a 
configuration with highest energy (the anti ground state). 
The round-trip time thus gives an estimate of the equi- 
libration time for the flat-histogram ensemble. For the 
Ising model the number of energy levels scales linearly 
with system size N. While the round-trip time of an 
unbiased random walker scales like r ~ TV 2 , it was re- 
cently shown that for various two-dimensional Ising mod- 
els the growth with the number of spins is significantly 
stronger for the biased flat-histogram random walker i& 
For the ferromagnetic and fully frustrated Ising models 
polynomial scaling, r ~ N 2A and r ~ TV 2 9 , was found, 
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FIG. 1: Convergence of the round-trip time r versus the mod- 
ification factor / in the Wang-Landau algorithm. Shown are 
the results for three randomly generated three-dimensional 
±J spin glass samples with TV = 3 3 = 27, 4 3 , 6 3 = 216 spins. 
The measurements are averaged over 500 independent runs 
for TV = 27, 64 and over 30 runs for TV = 216, respectively. 



FIG. 2: Left panels: Distribution of round-trip times r for 
5000 randomly generated spin glass samples of size TV = 3 3 = 
27 and TV = 5 3 = 125 respectively. Right panels: Distribution 
of the ratio of the number of first excited states to the num- 
ber of ground states g(Ei) / g(Eo) for the same system sizes. 
In all panels the solid lines indicate fits to fat-tailed Frechet 
extremal- value distributions. 



and exponential growth for the two-dimensional ± J spin 
glass 

For a given sample we find that the round-trip time 
measured during the iterations of the WL algorithm con- 
verges as the simulated ensemble approaches the flat- 
histogram ensemble. The convergence of round-trip 
times is illustrated for three randomly generated spin 
glass samples in Fig. ^ Since correct convergence to the 
round-trip times of the exact flat-histogram ensemble was 
shown for the two-dimensional ± J Ising spin glass^ we 
assume correct convergence for the 3D case as well and 
thereby justify that our results for the asymptotic round- 
trip times hold for any flat-histogram method. 



A. Sample-to-sample variation 

To study the sample dependence of the round-trip 
times we have analyzed 5000 randomly generated spin 
glass samples for TV = 3 3 , 4 3 , 5 3 , 6 3 and 1000 samples for 
TV = 8 3 = 512, respectively. To assure convergence of the 
measured round-trip times we restrict the measurement 
to the final step in the Wang-Landau algorithm. We find 
strong sample-to-sample variations over several orders of 
magnitude for fixed system size TV which is shown in the 
left panels of Fig. El (TV = 3 3 and TV = 5 3 ). For the 
spin glass of size TV = 8 3 = 512 the full distribution of 
round-trip times is shown in Fig. [21 The distribution cov- 
ers some 4 orders of magnitude and contains spin glass 
samples which were simulated between some minutes and 
about a month on a 500 MHz Pentium III CPU. 

To quantitatively analyze the extremal events in the 
tails of the distributions, we use extremal- value theoryiii 
The central limit theorem for extremal-value states that 
the extrema of random subsets of any distribution follow 
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FIG. 3: Log-log plot of the distribution of round-trip times 
t for the three-dimensional ± J spin glass with TV = 8 3 = 512 
spins. The round-trip times r were measured for the con- 
verged flat-histogram ensemble in the final step of the Wang- 
Landau algorithm. Data from 1000 randomly generated spin 
glass samples are shown. 



a generalized extremal-value distribution^ This distri- 
bution takes one of three characteristic forms: Frechet 
(algebraic decay, fat-tailed), Weibull (exponential decay) 
or Gumbel (faster than exponential decay, thin-tailed). 
Here we find that all measured round-trip times seem 
to follow a fat-tailed Frechet extremal-value distribution, 
similar to the two-dimensional ± J spin glass^ This im- 
plies that the central limit theorem applies even to the 
smallest possible subset - a single round-trip time. As 
a consequence, every single three-dimensional ±J spin 
glass sample constitutes an extremal event. The inte- 
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FIG. 4: Scaling of the shape parameter £ of the fitted Frechet 
distributions versus system size. The m-th moment of a fat- 
tailed Frechet distribution is well-defined only if £ < 1/m. 
The dashed lines indicate where variance and mean of the 
respective distributions become ill-defined. 



FIG. 5: Log-log plot of the running mean of round-trip times 
defined in Eq. © versus the number of samples used for 
averaging. Results for various system sizes are shown. For 
TV < 6 3 = 216 the mean is well-defined and the running mean 
seems to converge (dashed line). For TV = 8 3 = 512 the mean 
becomes ill-defined and the running mean diverges—. 



grated form of the Frechet distribution is given by 

H^At) = ex P {- (l + j ■ (3) 

The parameter /i indicates the location of the distribu- 
tion, that is the most probable round-trip time, and the 
parameter f3 defines the scale of the distribution, e.g. the 
height of the peak. The shape parameter £, which is posi- 
tive for fat-tailed distributions, describes the decay of the 
tail of the distribution. We have determined the three 
parameters /z, (3 and £ of the fitted Frechet distributions 
with a maximum likelihood estimator. The resulting fits 
are shown as solid lines in Fig. [21 

We now turn to the scaling of the shape parameter 
with system size TV which is shown in Fig. 0] With in- 
creasing system size the shape parameter monotonically 
increases. This strongly affects the fat tails of the distri- 
butions which in the limit r — > oo exhibit a power-law 
decay of the form 
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FIG. 6: Log-log plot of the running variance of the distribu- 
tion of round-trip times defined in Eq. ©. Results for various 
system sizes are shown. For TV < 4 3 = 64 the variance is well- 
defined and the running variance seems to converge (dashed 
line). For TV > 5 3 = 125 the variance becomes ill-defined and 
the running variance diverges—. 



A^m^t-u+vc). (4) 

From this asymptotic behavior we can see that the m-th 
moment of a fat-tailed Frechet distribution is well-defined 
only if £ < 1/m. The scaling in Fig. 0] suggests that 
for TV > 4 3 = 64 (TV > 8 3 = 512) the shape parameter 
becomes larger than 0.5 (1) and thus the variance (mean) 
of the distribution is no longer well-defined. 

To illustrate this unusual behavior we can calculate 
running moments of the distribution, e.g. by only con- 
sidering subsets of the first n round-trip times of all mea- 
sured round trip times {r} when calculating the moments 
of the distribution. The running mean of round-trip 



times is then defined by 

1 " 

Mean {r} (n) = - V n , (5) 

i=i 

and the running variance by 

1 n 2 

Var {r} (ra) = — — - (n - Mean {r} (n)) . (6) 
t=i 

The running mean and variance are calculated for a 
fixed random order of the set of round-trip times {r} = 
{n, T2, . . . , Tgooo}- Figs. and [S] show the running mean 
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FIG. 7: Correlation of the round-trip time r and ratios in the 
density of states for various system sizes. In the left panels 
the ratio Ri = g(Ei) / 'g(Eo) is shown. In the rights panels 
transitions to higher excited states are also included, R 2 = 
g(Ex)/g(E ) + g(E 2 )/g(E 1 ) + g(E 2 )/g(E ). Shown are data 
from 5000 randomly generated spin glass samples for N = 27, 
125, 216 and 1000 for N = 8 3 = 512, respectively. 



and variance for various system sizes. For small system 
sizes, N < 64, both mean and variance are well defined 
and the running mean and its error converges. For larger 
systems the variance becomes ill-defined as the shape pa- 
rameter £ turns larger than 0.5. Irregular "jumps" indi- 
cating rare events occur in the calculation of the running 
mean. For systems with N < 216 spins the running mean 
still converges, but the error of the running mean does 
not reduce even for large sample sets. For systems with 
more than TV = 8 3 = 512 spins the shape parameter be- 
comes larger than 1 and the mean thus ill-defined. This 
divergence of the mean round-trip time also becomes ap- 
parent in the calculation of the running mean where the 
irregular jumps occur so frequently that the mean round- 
trip time no longer converges^. 

This behavior becomes even more evident for the vari- 
ance, the second moment of the distribution, which ac- 
cording to the scaling of the shape parameter illustrated 
in Fig. 2| is no longer well-defined for systems larger than 
N > 4 3 = 64. As illustrated in Fig. the running 
variance diverges for larger systems^. Again, frequent 
"jumps" in the running variance indicate the occurrence 
of rare events in the fat tails of the distribution which 
dominate the calculation of the running variance. 



FIG. 8: Correlation of the computational effort, £ gC n.+CEA, 
of the genetic CEA (see text) and the ratio in the den- 
sity of states for N — 8 3 (1000 samples). Here the ratio 
R 2 = g(Ei)/g(E ) + g(E 2 )/g(E 1 ) + g(E 2 )/g(E ) is shown. 
The correlation to Ri = g(Ei)/g(Eo) looks similar. 



text of diffusivity measurements in section IIVI Here we 
study ratios in the density of states, such as the num- 
ber of first excited states to the number of ground states, 
g(Ei)/g(Eo). From the ground state a single-spin-flip up- 
date can connect at most Ng(Eo) states of energy E\ and 
thus the ratio g(Ei)/g(E ) gives a qualitative measure of 
the number of local minima, i.e. which are not reachable 
from a ground state via a single spin flip. In the right 
panels of Fig.^the distribution of this ratio g(E\) / g(Eo) 
is shown for 5000 spin glass samples of size N = 3 3 = 27 
and N = 5 3 = 125 respectively. Again we find that these 
distributions follow fat-tailed Frechet extremal- value dis- 
tributions. For these system sizes a strong correlation to 
the distribution of round-trip times is found which spans 
over several orders of magnitude as demonstrated in the 
left panels of Fig. For larger systems these correlations 
become less pronounced. However, if we consider addi- 
tional transitions, such as the transition from the second 
to first excited state, E2 — * E\ and E% — * Eq, and cal- 
culate the sum of the respective ratios in the density of 
states we can recover a correlation over several orders of 
magnitude as it is shown for systems with N = 6 3 = 216 
and N — 8 3 = 512 spins in the right panels of Fig. [7J 



C. Intrinsic correlations for the heuristic approach 



B. Intrinsic correlations for the WL algorithm 

In order to test whether the occurrence of the Frechet 
extremal-value distributions reflect an intrinsic property 
of the energy landscape three-dimensional spin glass we 
have analyzed the calculated density of states near the 
ground-state energy. The energy landscape near the 
ground state dominates to a large extent the measured 
round-trip times which is further discussed in the con- 



The strong correlation between intrinsic features of the 
energy landscape and the measured round-trip times for 
flat-histogram sampling naturally leads to the question 
whether the computational effort of other algorithms also 
comply with these intrinsic features. To this end, we com- 
pare the computational effort of the alternative heuristic 
approach with the density of states and the round-trip 
times measured in the WL algorithm for 1000 samples 
with N = 8 3 = 512 spins. 
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FIG. 9: Correlation of the computational effort, t gC n.+CEA, 
of the genetic CEA (see text) and the round-trip time t/N 2 
measured for multicanonical sampling in the final step of the 
WL algorithm for a system with TV = 8 3 spins (1000 samples). 

We define the computational effort of the genetic CEA 
as the running time of the algorithm which strongly de- 
pends on the parametrization of the underlying genetic 
algorithm, namely 

• Mi = initial size of population, i.e. how many con- 
figurations are optimized in parallel 

• rio = average number of offspring per configuration, 
i.e. how many iterations of the algorithm are run 

• ^min = number of CEA minimization steps per con- 
figuration and per iteration 

The running time then depends linearly on the product 
Mn n m i n . If the algorithm is run independently several 
times, one finds that not all runs result in ground states 
and we denote the fraction of runs which find the true 
ground state as /gs- In general one can observe that this 
fraction increases with the running time. 

We now define the computational effort for a 
given sample and parameter set {M{, n Q , rimin} a s 
i(Mi,n ,n min ) = Min n min / f GS (i.e. t = oo, if no 
ground state is found). For "simple" samples, many 
ground states are found for most parameter combina- 
tions, while "hard" samples need large sizes of popula- 
tions and/or many iterations and/or many minimization 
steps, thereby increasing the computational effort. For a 
given spin glass sample we perform multiple simulations 
with different combinations of parameters, and define the 
overall computational effort of a sample as the minimum 
over all parameter combinations considered (using a fixed 
value «r = 20): 

igon.+CEA = min {M { n n m - in / f GS ) . (7) 

(Mi,n ,n m i n ) 

In Fig. [S] the correlation of the computational effort 
of the genetic CEA algorithm is shown versus the ra- 
tio i?2 of the density of states as defined above for 1000 




FIG. 10: Scaling of the location (/x) and scale (/3) parame- 
ter of the Frechet distribution versus the linear system size 
L — 7V 1//3 . The parameters were fitted by a maximum- 
likelihood estimator. The solid lines indicate exponential fits 
of the respective data points. 



samples with N = 8 3 = 512 spins. There is only a 
weak correlation indicating that the heuristic algorithm is 
less sensitive to the energy-landscape close to the ground 
state than the WL algorithm. A direct comparison be- 
tween the genetic CEA and the WL algorithm is shown 
in Fig. [§] We find that the computational effort for the 
heuristic algorithm spreads over only 2 orders of magni- 
tude in comparison to some 4 orders of magnitude for 
the WL algorithm. The correlation between the two al- 
gorithms is weak, and less pronounced for those sam- 
ples which are especially hard to equilibrate using multi- 
canonical sampling. 



D. Scaling of typical round-trip times 

Although the mean round-trip times are no longer well- 
defined for larger systems, the location and scale param- 
eter of the Frechet distribution stay well-defined and can 
be used to further characterize the scaling of these fat- 
tailed distributions with system size. Fig. 1101 shows that 
these parameters exponentially diverge with linear sys- 
tem size like 

H oc exp(Z/(1.34±0.03)) , 

P oc exp(Z/(1.00±0.03)) . (8) 

The study of the equilibrium behavior of larger system 
size is thus not only limited by the occurrence of rare 
events, but also by the exponential growth of the round- 
trip times in the "bulk" of the distribution which ren- 
ders a comprehensive study of the sample-to-sample vari- 
ations impossible. 
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FIG. 11: Local diffusivity D(E,to) in energy of a flat- 
histogram random walker for the three-dimensional Ising fer- 
romagnet. Data for various diffusion times are shown for a 
system with N — 4 3 = 64 spins. 



IV. DIFFUSIVITY MEASUREMENTS 

We have seen that both the asymptotic and the dy- 
namic performance are limited by the rough energy land- 
scape close to the ground-state. To study these limita- 
tions in more detail we measure the local diffusivity of 
the random walker in energy. Recently, it was shown 
that the simulated statistical ensemble can be optimized 
by a feedback loop which reweights the ensemble based 
on preceeding measurements of the local diffusivity^ Al- 
though we do not follow up on this idea in the present 
study, we can measure the diffusivity to reveal the "bot- 
tlenecks" of the biased random walk in energy as local 
minima in the diffusivity. Here we use a time-dependent 
definition of the diffusivity, D(E,tD), in energy space 



D(E, t D ) = {{E{t) - E(t + t D )) 2 ^ /t D , 



(9) 



where to is the diffusion time. The relevant time scale for 
the diffusion time is set by the round-trip time in energy. 

As a first example we study the local diffusivity of 
the three-dimensional ferromagnetic Ising model which 
undergoes a second-order phase transition from a mag- 
netically ordered to a disordered phase at finite energy 
E C /3N ^ -0.43,22 The measured diffusivity D(E,t D ) of 
the flat-histogram random walker is shown in Fig.^] We 
find that the diffusivity is not constant as it is expected 
for an unbiased random walk, but there is a broad mini- 
mum below the critical energy. In this energy region the 
random walker is slowed down due to the slow dynamics 
of domain walls which separate droplets of magnetically 
ordered phases. 

Next we turn to the three-dimensional spin glass. Mea- 
surements of the diffusivity of the flat-histogram random 
walker for two randomly generated samples are shown 
in Fig. 1121 For both samples we find a minimum of the 
diffusivity at the ground-state energy. We further note 
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FIG. 12: Local diffusivity D(E,to) in energy of a flat- 
histogram random walker for two three-dimensional ± J Ising 
spin glass samples. Data for various diffusion times are shown 
for systems with N = 4 3 = 64 spins. 



that with increasing round-trip times the minimum in 
the diffusivity becomes more pronounced. These diffusiv- 
ity measurements further underline that the bottleneck 
of the flat-histogram random walker are at the ground- 
state energy. The suppressed diffusivity in this energy re- 
gion gives rise to an entropic barrier which results in long 
round-trip times and aggravates equilibration of the ran- 
dom walker in the low energy phase. For larger systems 
the suppressed local diffusivity of the flat-histogram ran- 
dom walker renders a comprehensive study of the glassy 
phase impossible as it becomes computationally too ex- 
pensive to equilibrate a large number of samples needed 
to calculate statistical averages. The recent development 
of a systematic means to optimize the simulated statisti- 
cal ensemble^ holds promise to overcome some of these 
limitations. 



V. CONCLUSIONS 

To summarize, we have studied the performance of 
the Wang-Landau algorithm for the three-dimensional 
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±J Ising spin glass. The asymptotic performance 
- which corresponds to the performance of any flat- 
histogram method sampling a multicanonical ensemble 
such as the multicanonical method, 4 simulated and par- 
allel tempering^ broad histograms^ and transition ma- 
trix Monte Carloi — is found to be dominated by strong 
sample-to-sample variations. The measured round-trip 
times follow fat-tailed Frechet extremal-value distribu- 
tions. The typical round-trip times in the bulk of these 
distributions obey exponential scaling, as expected for 
this NP-hard problem. The fat tails of the distributions 
dominate the calculation of statistical averages. A care- 
ful statistical analysis is needed which goes beyond the 
calculation of the moments of a finite sample distribution 
as it was done in previous studies. The intrinsic character 
of the Frechet extremal distributions becomes evident in 
strong correlations over several orders of magnitude be- 
tween the round-trip time and the behavior of the density 
of states near ground-state energy. The origin of the ex- 
tremal character of every single spin glass sample remains 
an open question which deservers further investigations. 

Our measurements of the dynamic performance and 
comparison with ground states obtained by genetic CEA 
showed that for samples with up to N — 8 3 spins the 
Wang-Landau algorithm always finds the correct ground- 
state energies. For sizes 10 3 , 12 3 , one can find true 
ground states, if one restricts the random walk to a small 
energy bin around the exact ground-state energy calcu- 
lated by the heuristic CEA approach. For samples with 
N = 1 4 3 spins we identified samples for which the Wang- 



Landau algorithm does not find the true ground-state 
energy within reasonable simulation times (« 10 7 MC 
sweeps) and does not converge towards the multicanoni- 
cal ensemble. 

The cntropic barrier near the zero-temperature ground 
state can be revealed as a pronounced minimum in the 
local diffusivity of the flat-histogram random walker in 
energy. The recent development of optimized statisti- 
cal ensembles based on feedback of the local diffusivity 
holds promise to overcome the observed slow down of 
the flat-histogram random walker and thereby enhancing 
equilibration in these systems^ 
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